#if !defined( SPARSE_MATRIX_H)
#define  SPARSE_MATRIX_H
#include <algorithm>
#include <memory>
#include <iostream>
#include <vector>
#include <string>
#include <cassert>
#include <functional>
#include "stlhelper.h"
#include <bits/stl_function.h>
#include <boost/scoped_array.hpp>
#include <boost/numeric/ublas/matrix_sparse.hpp>


#include "matrix_traits.h"
#include "matrix_common.h"
#include "matrix.h"

class t_progress_indicator;
namespace bcpp_matrix
{

//________________________________________________________________________________________
//
//  Sparse matrix implementation
//
//      which contains
//          row/column names
//          matrix names
//
//
//________________________________________________________________________________________

template <typename T> 
struct t_sparse_matrix
{
public:
	typedef typename boost::numeric::ublas::compressed_matrix<T> sparse_matrix_type;
	//compressed_matrix<double> m (11, 11, 3);
	
	sparse_matrix_type			matrix;
	std::string					matrix_name;
	std::vector<std::string>	names;


	// symmetrical matrix
	unsigned					size() const {return matrix.size1();}
	//t_sparse_matrix<T>(t_sparse_matrix& 	matrix_):
    //                                      matrix(matrix_)
	t_sparse_matrix<T>()					{}
	t_sparse_matrix<T>(unsigned sz)			{allocate(sz);}
	~t_sparse_matrix(){}
	t_sparse_matrix<T>(const t_matrix<T>& other);
	t_sparse_matrix<T>(const t_dense_matrix<T>& other, float null_value);
	
	// don't assume asymmetrical matrix
	unsigned					get_count_filled() const
										{ return matrix.nnz();}
	unsigned					get_count_pairs() const
										{ return size () * size();}
	unsigned					get_count_missing() const
										{ return get_count_pairs() - get_count_filled();}
	bool						is_missing(unsigned i,  unsigned j) const
										{return !matrix.find_element(i, j);	}
	void						erase_element(unsigned i,  unsigned j)
										{matrix.erase_element(i, j);}

	void allocate(unsigned size)
	{
		matrix.resize(size, size, false);
		names.assign(size, std::string());
	}
	void do_count_filled();
	void swap(t_sparse_matrix& other, const std::vector<std::string>& names_)
	{
		names = names_;
		matrix.swap(other);
	}


	void swap(t_sparse_matrix<T>& other)
	{
		matrix.swap(other.matrix);
		names.swap(other.names);
		std::swap(matrix_name, other.matrix_name);
	}
	void get_single_linkage_clusters(std::vector<t_sparse_matrix<T> >& data) const;
	void reduce_to_subset(const std::vector<unsigned>& order_indices);
	void subset(const std::vector<unsigned>& indices, t_sparse_matrix<T>& new_data) const;
	void sort_by_name_alpha();
	
};

template <typename T> 
bool operator == (const t_sparse_matrix<T>& a, const t_sparse_matrix<T>& b)
{
	// DEBUGG
	if (a.matrix_name != b.matrix_name)
		std::cerr << "\nmatrix name!=" << a.matrix_name << "!=" << b.matrix_name << "\n";
	if (a.names.size() != b.names.size())
		std::cerr << "\nmatrix names!=" << a.names.size() << "!=" << b.names.size() << "\n";
	if (a.matrix.nnz() != b.matrix.nnz())
		std::cerr << "\nmatrix sizes!=\n";
	if (a.names != b.names)
		std::cerr << "\nmatrix names1!=\n";

	// compare name
	if (a.matrix_name != b.matrix_name)
		return false;

	// compare size
	if (a.names.size() != b.names.size())
		return false;
	if (a.matrix.size1() != b.matrix.size1())
		return false;
    //if (a.matrix.size2() != b.matrix.size2())
    //  return false;

	// compare filled size
	if (a.matrix.nnz() != b.matrix.nnz())
		return false;

	// compare column names
	if (a.names != b.names)
		return false;


	// compare matrix. All filled positions must be equal and filled.
    // since already checked filled sizes are equal, matrices must
    // aso be equal
	typedef typename t_sparse_matrix<T>::sparse_matrix_type::const_iterator1 t_i1;
	typedef typename t_sparse_matrix<T>::sparse_matrix_type::const_iterator2 t_i2;
	typedef const T* p_T;
	
	for (t_i1 i1 = a.matrix.begin1(), end1 = a.matrix.end1(); i1 != end1; ++i1)
	{
		for (t_i2 i2 = i1.begin(); i2 != i1.end(); ++i2)
		{
			unsigned ii = i2.index1();
			unsigned jj = i2.index2();
			p_T p_b_value = b.matrix.find_element(ii, jj);
			if (!p_b_value)
				std::cerr << "\nMissing "<< ii << ", " << jj << "\n";
			if (p_b_value && *i2 != *p_b_value)
				std::cerr << "\nDiff at "<< ii << ", " << jj << ":-"
					<< *i2 << "!=" << *p_b_value << "\n";
			if (!p_b_value || *i2 != *p_b_value)
				return false;
		}
	}
	return true;
}
template <typename T> 
inline void swap(t_sparse_matrix<T>& lhs, t_sparse_matrix<T> rhs) {lhs.swap(rhs);}


typedef t_sparse_matrix<float> t_sparse_matrix_f;


//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

//	get part of matrix 

//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

template<typename T>
typename t_sparse_matrix<T>::sparse_matrix_type&	
get_sub_matrix(	const std::vector<unsigned>&				indices,
				const typename t_sparse_matrix<T>::sparse_matrix_type&	original_matrix,
				typename t_sparse_matrix<T>::sparse_matrix_type&		new_matrix)
{
	if (!indices.size())
		return new_matrix;
	typedef std::vector<unsigned>::const_iterator t_uiter;
	new_matrix.resize(indices.size(), indices.size(), false);

	typedef const T* p_T;
	for (unsigned i = 0; i < indices.size(); ++i)
	{
		for (unsigned j = 0; j < indices.size(); ++j)
		{
			assert(indices[i] < original_matrix.size1());
			assert(indices[j] < original_matrix.size1());
			p_T p_orig_element = original_matrix.find_element(indices[i], indices[j]);
			if (p_orig_element)
				new_matrix(i, j) = *p_orig_element;
		}
	}
	return new_matrix;
}




//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

//	single linkage clusters

//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888


template <class T>
void t_sparse_matrix<T>::get_single_linkage_clusters(std::vector<t_sparse_matrix<T> >& data) const
{

	//
	//	separate matrix using single linkage clustering
	//
	std::vector<std::vector<unsigned> > clusters;
	private_details::get_single_linkage_clusters_indices(*this, clusters);
//	if (clusters.size() == 1)
//	{
//		data.push_back(t_sparse_matrix<T>());
//		data.back() = *this;
//		return;
//	}

	//
	//	for each cluster
	//
	for (unsigned i = 0; i < clusters.size(); ++i)
	{
		data.push_back(t_sparse_matrix<T>());

		// divide up data
		get_sub_matrix<T>(clusters[i], matrix, data.back().matrix);

		// divide up names
		private_details::extract_clustered_names(names , clusters[i], data.back().names);
		
		// data' name is my name
		data.back().matrix_name = matrix_name;
	}
}




template <typename T> 
void t_sparse_matrix<T>::reduce_to_subset(const std::vector<unsigned>& order_indices)
{
	t_sparse_matrix<T> new_data;
	subset(order_indices, new_data);
	swap(new_data);
}

template <typename T> 
void t_sparse_matrix<T>::subset(const std::vector<unsigned>& indices, t_sparse_matrix<T>& new_data) const
{
	// copy values to new data
	new_data.matrix_name= matrix_name;

	// get data subset
	get_sub_matrix<T>(indices, matrix, new_data.matrix);

	// get names subset
	for (unsigned i = 0; i < indices.size(); ++i)
		new_data.names.push_back(names[indices[i]]);
}



template <typename T> 
void t_sparse_matrix<T>::sort_by_name_alpha()
{
	std::vector<unsigned> indices;
	sort_indices(names, indices);
	reduce_to_subset(indices);
}

template <typename T> 
t_sparse_matrix<T>::t_sparse_matrix(const t_matrix<T>& other)
	:
	matrix(other.matrix_name),
	names(other.names)
{
	matrix.allocate(other.size());
	unsigned size = other.size();
	for (unsigned i = 0; i < size; ++i)
		for (unsigned j = 0; j < size; ++j)
			if (!other.matrix.is_missing(i, j))
				matrix(i, j) = other.matrix[i][j];
}
template <typename T> 
t_sparse_matrix<T>::t_sparse_matrix(const t_dense_matrix<T>& other, float null_value)
{
	allocate(other.size());
	unsigned size = other.size();
	for (unsigned i = 0; i < size; ++i)
		for (unsigned j = 0; j < size; ++j)
			if (other[i][j] != null_value)
				matrix(i, j) = other[i][j];
}

};	// bcpp_matrix

	

#endif	//SPARSE_MATRIX_H

